#include <iostream>
#include <fstream>
#include <ctime>
#include <cstdlib>
#include <cmath>
#include "anyoption.h"
#include "sbml/SBMLTypes.h"
#include "../libBioModel/Network.h"
#include "../libBioModel/Stoichiometry.h"
#include "../libBioModel/State.h"
#include "../libBioModel/InitialConds.h"

using namespace std;


/*
 * An implementation of Gillespie's Stochastic Simulation Algorith to solve 
 * the Chemical Master Equation based on Monte Carlo simulation.
 *
 * This program implements the direct method of SSA.
 *
 * Sep 22, 2009
 * Youfang Cao
 *
 */

int main(int argc, char* argv[])
{
	/* Variable definitions */
	char*		sbmlfile = 0;        	// SBML filename.
	char*		initfile = 0;        	// Inital condition filename.
	unsigned int	n_timesteps = 100000;    		// The number of time steps.
	double		totaltime = 1000;            	// Total time for simulation.
	int		verbose = 0;
	char*		outputfile = 0;


	// Iterative variables
	unsigned int i=0, j=0, k=0; // some index variables.
	int f = 0;
	unsigned int r = 0;
	double t = 0;
	double dt = 0;

	// Variable definitions for Gillespie Algorithm
	InitialConds X0; // The initial state.
	State* X = 0; // The state.
	Network nw;
	unsigned int n_reactions = 0;
	double tau1 = 0;
	double tau2 = 0;
	double* rates = 0;
	Reaction* re;
	KineticLaw* klaw;
	double a0 = 0;
	double cum = 0;
	double cri = 0;

        /* CREATE AN OBJECT TO HANDLE COMMAND LINE */
        AnyOption *opt = new AnyOption();

        /* SET THE USAGE/HELP   */
        opt->addUsage( "" );
        opt->addUsage( "Usage: " );
        opt->addUsage( "" );
        opt->addUsage( " -h  --help  		Prints this help " );
        opt->addUsage( " -m  --sbml  		Specify the SBML file" );
        opt->addUsage( " -i  --init  		Specify the initial condition" );
        opt->addUsage( " -S  --nsteps 100000    Number of time steps for each simulation " );
        opt->addUsage( " -t  --ttime 1000  	Total time for simulation " );
        opt->addUsage( " -o  --ofile 		Output file " );
        opt->addUsage( " -v  --verbose 0  	Switch for more output " );
        opt->addUsage( "" );

        /* SET THE OPTION STRINGS/CHARACTERS */
        opt->setFlag(  "help", 'h' );
        opt->setOption(  "sbml", 'm' );
        opt->setOption(  "init", 'i' );
        opt->setOption(  "nsteps", 'S' );
        opt->setOption(  "ttime", 't' );
        opt->setOption(  "ofile", 'o' );
        opt->setOption(  "verbose", 'v' );

	/* go through the command line and get the options  */
        opt->processCommandArgs( argc, argv );

        /* GET THE VALUES */
        if( opt->getFlag( "help" ) || opt->getFlag( 'h' ) ) 
	{
                opt->printUsage();
		delete opt;
		return 1;
	}
	if( opt->getValue( 'm' ) != NULL  || opt->getValue( "sbml" ) != NULL  )
		sbmlfile = opt->getValue( 'm' );
	else
	{
		printf ("Please specify the SBML filename.\n\n");
		opt->printUsage();
		return 1;
	}
	if( opt->getValue( 'i' ) != NULL  || opt->getValue( "init" ) != NULL  )
		initfile = opt->getValue( 'i' );
	else
	{
		printf ("Please specify the initial condition.\n\n");
		opt->printUsage();
		return 1;
	}
	if( opt->getValue( 'S' ) != NULL  || opt->getValue( "nsteps" ) != NULL  )
		n_timesteps = atoi(opt->getValue( 'S' ));
	if( opt->getValue( 't' ) != NULL  || opt->getValue( "ttime" ) != NULL  )
		totaltime = atof(opt->getValue( 't' ));
	if( opt->getValue( 'v' ) != NULL  || opt->getValue( "verbose" ) != NULL  )
		verbose = atoi(opt->getValue( 'v' ));
	if( opt->getValue( 'o' ) != NULL  || opt->getValue( "ofile" ) != NULL  )
		outputfile = opt->getValue( 'o' );
	else
	{
		printf ("Please specify the output filename.\n\n");
		opt->printUsage();
		return 1;
	}
	cout << endl;

	/* Initialize variables. */

	/* READ SBML MODEL */
	nw.d = readSBML(sbmlfile);
	nw.m = nw.d->getModel();

	//errors = nw.d->getNumErrors();
	
	/* read stoichiometic matrix */
	nw.stoichiometry = new Stoichiometry();
	f = nw.stoichiometry->genStoichiometricMatrix(nw.m);
	if (f == -1)
	{
		printf ("Exit caused by error.\n");
		return -1;
	}
	if (verbose > 0)
		nw.stoichiometry->printStoichiometricMatrix();
	
	//modelid = nw.m->getId();
	nw.n_compartments = nw.m->getNumCompartments();
	nw.n_species = nw.m->getNumSpecies();
	nw.n_rules = nw.m->getNumRules();
	n_reactions = nw.stoichiometry->getNumReactions();
	nw.n_reactions = nw.m->getNumReactions();

	/* read reaction IDs */
	nw.ReactionsID = new string[n_reactions];
	j = 0;
	for (i = 0; i < nw.n_reactions; i++)
	{ nw.ReactionsID[j] = nw.m->getReaction(i)->getId();
		//cout << "heloo. " << n_reactions << " " << i << " " << j << " " << nw.ReactionsID[j] << endl;
		j++;
		if (nw.m->getReaction(i)->getReversible())
		{
			nw.ReactionsID[j] = nw.m->getReaction(i)->getId();
			j++;
		}
	}

	//cout << "heloo. " << endl;
	/* read compartment IDs */
	nw.CompartmentsID = new string[nw.n_compartments];
	for (i = 0; i < nw.n_compartments; i++)
	{
		nw.CompartmentsID[i] = nw.m->getCompartment(i)->getId();
	}

	/* read species IDs */
	nw.SpeciesID = new string[nw.n_species];
	if (verbose > 0)
	{
		printf ("Species:\n");
		printf ("id\tinitAmt\tinitCon\tconst\tboundcond\n");
	}
	for(i = 0; i < nw.n_species; i++)
	{
		/* Store id data in the array. */
		nw.SpeciesID[i] = nw.m->getSpecies(i)->getId();
	}
	
	/*
	 * Get rules ID list.
	 */
	if(verbose > 0)
	{
		printf ("Rules:\n");
	}
	nw.RulesID = new string[nw.n_rules];
	nw.n_assignmentRules = 0;
	for(i = 0; i < nw.n_rules; i++)
	{
		if (SBase_getTypeCode((SBase_t*)nw.m->getRule(i)) != SBML_ASSIGNMENT_RULE)
			continue;
			
		/* Store id data in the array. */
		nw.RulesID[i] = ((AssignmentRule*)nw.m->getRule(i))->getVariable();
		
		nw.n_assignmentRules ++;
	}

	/*
	 * Get parameters ID list.
	 */
	if(verbose > 0)
	{
		printf ("Global Parameters:\n");
	}
	nw.n_parameters = nw.m->getNumParameters();
	nw.ParametersID = new string[nw.n_parameters];
	for (i = 0; i < nw.n_parameters; i++)
	{
		nw.ParametersID[i] = nw.m->getParameter(i)->getId();
		if (verbose > 0)
		{
			printf ("%s\t%f\n", nw.ParametersID[i].data(), nw.m->getParameter(i)->getValue());
		}
	}

	/* Local parameter IDs */
	nw.LocalParamIDs = new string[20];

	/* Initialize rates */
	n_reactions = nw.stoichiometry->getNumReactions();
	rates = new double[n_reactions];
	for (i=0; i<n_reactions; i++)
	{
		rates[i] = 0;
	}

	/* READ INITIAL CONDITION */
	if (X0.readInitialConds(initfile) == -1)
	{
		printf ("ERROR:  Read initial condition file error.\n");
		return -1;
	}



	/* START SSA ITERATIONS */
	ofstream ftrajectory;
	ftrajectory.open(outputfile);
	for (i = 0; i < nw.n_species; i++)
	{
		ftrajectory << nw.SpeciesID[i] << " ";
	}
	ftrajectory << endl;

	srand(time(NULL));
	i = 0;
	t = 0;
	X = new State();
	X->CopyfromState(X0.getInitialState(0));
	while (i < n_timesteps) // && t < totaltime)
	{
		//cout << "Run: " << i+1;
		X->PrintState(ftrajectory);

		/* Calculate transition rates for all reactions in current state X */
		a0 = 0;
		for (k=0; k<n_reactions; k++)
		{
			re = nw.m->getReaction(abs(nw.stoichiometry->getReactionAttribute(k)) - 1);
			klaw = re->getKineticLaw();
			rates[k] = nw.getReactionRate(klaw, k, X) * (double)(nw.stoichiometry->canReact(k, X));
			a0 = a0 + rates[k];
			//cout << rates[k] << " ";
		}
		//cout << endl;

		/* Find the tau1 and tau2 */
		tau1 = (double)rand()/(double)RAND_MAX;
		tau2 = (double)rand()/(double)RAND_MAX;

		/* Now, I got the delta t */
		dt = log(1/tau1)/a0;

		//cout << dt << " " << tau1 << " " << a0 << " " << tau2 << " " << a0*tau2 << endl;

		cum = 0;
		cri = tau2*a0;
		for (k=0; k<n_reactions; k++)
		{
			cum += rates[k];
			if (cum >= cri) { break; }
		}
		/* Now, I got which reaction to go */
		r = k;

		/* Update the system state X and time t */
		t = t + dt;
		nw.stoichiometry->updateStateWithReaction(r, X);


		/* Go ahead to next run */
		i ++;
	}



}


/*
 */




